################################################################################
# PACKAGE AND DATA LOADING
################################################################################
devtools::install_github("jpmam1/sensemakr")

pacman::p_load(dplyr,ggplot2,ggcorrplot,fixest,modelsummary,sensemakr,AER,ivDiag,did,
               stargazer,tidyr,scales)

df <- readRDS('master_df.rds')
year_counts <- readRDS("walls_year.rds")

################################################################################
# FIGURE 1
################################################################################

png("Figure 1.png",4000,3500,res=500)
ggplot(df, aes(year, torture_vdem)) +
  geom_smooth(method = "loess", colour="violet", size=1) + ylab("Freedom from Torture Index\n(V-Dem, inversed)") +
  xlab("Year") + theme_bw() +
  labs(caption = "Higher values in the Freedom of Torture Index represent\nmore practices of torture by state officials.")
dev.off()

################################################################################
# FIGURE 2
################################################################################

# Set theme
theme_fig2 <- function(base_size = 11, base_family = "serif") {
  theme_minimal(base_size = base_size, base_family = base_family) +
    theme(text = element_text(color = "black"),
          plot.title = element_text(face = "bold",size = base_size + 1,hjust = 0.5,margin = margin(b = 10)),
          plot.subtitle = element_text(hjust = 0.5,color = "gray40", margin = margin(b = 15)),
          axis.title = element_text(face = "bold"),
          axis.text = element_text(color = "black"),
          axis.line = element_line(color = "black", linewidth = 0.5),
          axis.ticks = element_line(color = "black", linewidth = 0.5),
          panel.grid.major = element_line(color = "gray90", linewidth = 0.3),
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "white", color = NA),
          plot.background = element_rect(fill = "white", color = NA),
          legend.position = "bottom",
          legend.box = "horizontal",
          legend.title = element_text(face = "bold"),
          legend.text = element_text(size = base_size - 1),
          legend.margin = margin(t = 10),
          plot.margin = margin(20, 20, 20, 20),
          strip.text = element_text(face = "bold", size = base_size),
          strip.background = element_rect(fill = "gray95", color = NA)
    )
}

# Set default theme
theme_set(theme_fig2(base_size = 11, base_family = "Times"))

df_figure2 <- df %>% group_by(year) %>%
  summarise(torture_bio=mean(amelia_fit,na.rm=TRUE))
df_figure2 <- merge(df_figure2,year_counts,by="year",all.x=TRUE)

range_torture <- range(df_figure2$torture_bio, na.rm = TRUE)
range_walls <- range(df_figure2$walls, na.rm = TRUE)
scale_factor <- diff(range_torture) / diff(range_walls)
shift <- range_torture[1] - range_walls[1] * scale_factor
df_plot <- df_figure2 %>%
  mutate(walls_scaled = walls * scale_factor + shift,year_factor = factor(year))

p <- ggplot(df_plot) +
  geom_line(aes(x = year, y = torture_bio, linetype = "Level of Torture by BIOs"),
    linewidth = 0.8,color = "black") +
  geom_point(aes(x = year, y = torture_bio, shape = "Level of Torture by BIOs"),size = 2.5,
    color = "black",fill = "white") +
  geom_line(aes(x = year, y = walls_scaled, linetype = "Border Walls"),linewidth = 0.8,
            color = "gray40") +
  geom_point(aes(x = year, y = walls_scaled, shape = "Border Walls"),size = 2.5,color = "gray40",
    fill = "white") +
  scale_x_continuous(name = "Year",breaks = pretty_breaks(),expand = expansion(mult = 0.05)) +
  scale_y_continuous(name = "Level of Torture by BIOs",
                     sec.axis = sec_axis(~ (. - shift) / scale_factor,name = "Border Walls"),
    breaks = pretty_breaks(),expand = expansion(mult = 0.05)) +
  scale_shape_manual(name = NULL,values = c("Level of Torture by BIOs" = 21, "Border Walls" = 24)) +
  scale_linetype_manual(name = NULL,values = c("Level of Torture by BIOs" = "solid", "Border Walls" = "dashed")) +
  labs(x = "Year")

png("Figure 2.png",3000,2000,res=500)
p
dev.off()

################################################################################
# COEFFICIENT MAP DEFINITIONS
################################################################################

cm_h1 <- c(
  "torture_vdem"       = "Torture (V-Dem)",
  "tort"               = "Torture (CIRI)",
  "physint"            = "Human rights (Fariss)",
  "rv"                 = "Random Variable",
  "Police"             = "Police",
  "Prison"             = "Prison",
  "Military"           = "Military",
  "Intelligence"       = "Intelligence",
  "Paramilitary"       = "Paramilitary",
  "liberal_democracy"  = "Liberal democracy",
  "ideology"           = "Ideology: Right",
  "hihost"             = "Inter-state Conflict",
  "civil_war"          = "Civil War",
  "log(gdp_pc)"        = "GDP pc (ln)",
  "log(population)"    = "Population (ln)"
)

cm_main <- c(
  "wall_const"                         = "Wall",
  "fit_wall_const"                     = "Wall",
  "wall_const3"                        = "Wall (3 years)",
  "log(nonfinancial + 8.04001765049333)" = "Investment",
  "structure_tot"                      = "Border Infrastructure",
  "visa"                               = "Visa Policy",
  "visa_extended"                      = "Visa Policy",
  "v2csreprss"                         = "CSO Repression",
  "wall_const:v2csreprss"              = "Wall x CSO Repression",
  "wall_const:visa_extended"           = "Wall x Visa Policy",
  "mean_anxiety"                       = "Anxiety",
  "mean_anger"                         = "Anger",
  "wall_const:mean_anxiety"            = "Wall x Anxiety",
  "wall_const:mean_anger"              = "Wall x Anger",
  "wall_const3:mean_anxiety"           = "Wall (3 years) x Anxiety",
  "wall_const3:mean_anger"             = "Wall (3 years) x Anger",
  "structure_tot:mean_anxiety"         = "Border Infrastructure x Anxiety",
  "structure_tot:mean_anger"           = "Border Infrastructure x Anger",
  "visa_extended:mean_anxiety"         = "Visa Policy x Anxiety",
  "visa_extended:mean_anger"           = "Visa Policy x Anger",
  "log(migration_intpl + 1)"           = "Migrants",
  "wall_const:log(migration_intpl + 1)"= "Wall x Migrants",
  "log(refugees + 1)"                  = "Refugees",
  "wall_const:log(refugees + 1)"       = "Wall x Refugees",
  "log(asylum_seekers + 1)"            = "Asylum seekers",
  "civil_war"                          = "Civil War",
  "log(fatalities + 1)"                = "Violence (neighbor)",
  "ideology"                           = "Ideology: Right",
  "torture_vdem"                       = "Torture",
  "liberal_democracy"                  = "Liberal democracy",
  "wall_reporting"                     = "Wall reporting",
  "log(wall_reporting + 1)"            = "Wall reporting",
  "log(population)"                    = "Population (log)",
  "log(nonfinancial + invt_min_1)"     = "Investment"
)

################################################################################
# HYPOTHESIS 1 — BIOs TORTURE VS GENERAL TORTURE PATTERNS
################################################################################

# --- Table 1 models ---
h1_1  <- feols(BIOs ~ torture_vdem | ccode + year, data = df, cluster = "ccode")
h1_2  <- feols(BIOs ~ tort | ccode + year, data = df, cluster = "ccode")
h1_3  <- feols(BIOs ~ physint | ccode + year, data = df, cluster = "ccode")
h1_4  <- feols(BIOs ~ torture_vdem + liberal_democracy + log(population) + log(gdp_pc) + hihost + civil_war + ideology | ccode + year, data = df, cluster = "ccode")
h1_5  <- feols(BIOs ~ physint + torture_vdem + liberal_democracy + log(population) + log(gdp_pc) + hihost + civil_war + ideology | ccode + year, data = df, cluster = "ccode")
h1_6  <- feols(BIOs ~ Police + Prison + Military + Intelligence + Paramilitary | ccode + year, data = df, cluster = "ccode")
h1_7  <- feols(BIOs ~ Police + Prison + Military + Intelligence + Paramilitary + liberal_democracy + log(population) + log(gdp_pc) + hihost + civil_war + ideology | ccode + year, data = df, cluster = "ccode")
h1_16 <- feols(BIOs ~ torture_vdem + liberal_democracy + log(population) + log(gdp_pc) + hihost + civil_war + ideology | year, data = df, cluster = "ccode")
h1_18 <- feols(BIOs ~ Police + Prison + Military + Intelligence + Paramilitary | year, data = df, cluster = "ccode")
h1_19 <- feols(BIOs ~ Police + Prison + Military + Intelligence + Paramilitary + liberal_democracy + log(population) + log(gdp_pc) + hihost + civil_war + ideology | year, data = df, cluster = "ccode")

# --- Table A.2 models (bivariate by agent type) ---
h1_8  <- feols(BIOs ~ Police | ccode + year, data = df, cluster = "ccode")
h1_9  <- feols(BIOs ~ Prison | ccode + year, data = df, cluster = "ccode")
h1_10 <- feols(BIOs ~ Military | ccode + year, data = df, cluster = "ccode")
h1_11 <- feols(BIOs ~ Intelligence | ccode + year, data = df, cluster = "ccode")
h1_12 <- feols(BIOs ~ Paramilitary | ccode + year, data = df, cluster = "ccode")

# --- Table A.3 models (neighbor covariates) ---
h1_n1 <- feols(BIOs ~ torture_vdem + liberal_democracy + log(population_neigh_total) + log(gdp_pc_neigh_mean) + hihost + civil_war_neigh + ideology | ccode + year, data = df, cluster = "ccode")
h1_n2 <- feols(BIOs ~ physint + torture_vdem + liberal_democracy + log(population_neigh_total) + log(gdp_pc_neigh_mean) + hihost + civil_war_neigh + ideology | ccode + year, data = df, cluster = "ccode")
h1_n3 <- feols(BIOs ~ Police + Prison + Military + Intelligence + Paramilitary + liberal_democracy + log(population_neigh_total) + log(gdp_pc_neigh_mean) + hihost + civil_war_neigh + ideology | ccode + year, data = df, cluster = "ccode")
h1_n4 <- feols(BIOs ~ torture_vdem + liberal_democracy + log(population_neigh_mean) + log(gdp_pc_neigh_mean) + hihost + civil_war_neigh + ideology | ccode + year, data = df, cluster = "ccode")
h1_n5 <- feols(BIOs ~ physint + torture_vdem + liberal_democracy + log(population_neigh_mean) + log(gdp_pc_neigh_mean) + hihost + civil_war_neigh + ideology | ccode + year, data = df, cluster = "ccode")
h1_n6 <- feols(BIOs ~ Police + Prison + Military + Intelligence + Paramilitary + liberal_democracy + log(population_neigh_mean) + log(gdp_pc_neigh_mean) + hihost + civil_war_neigh + ideology | ccode + year, data = df, cluster = "ccode")

# --- Random variable placebo (for Table A.1) ---
h1_1_rv <- feols(BIOs ~ rv | ccode + year, data = df, cluster = "ccode")
h1_4_rv <- feols(BIOs ~ rv + liberal_democracy + log(population) + log(gdp_pc) + hihost + civil_war + ideology | ccode + year, data = df, cluster = "ccode")
h1_5_rv <- feols(BIOs ~ physint + rv + liberal_democracy + log(population) + log(gdp_pc) + hihost + civil_war + ideology | ccode + year, data = df, cluster = "ccode")

# --- Table A.1: Goodness-of-fit comparison ---
adj <- function(m) 1 - (((m$nobs - 1) * (1 - m$sq.cor)) / (m$nobs - m$nparams - 1))

test_h1 <- data.frame(
  Model = c("Adj. R2 (Table 1)", "Adj. R2 (Random variable)"),
  `(1)` = c(adj(h1_1),    adj(h1_1_rv)),
  `(2)` = c(adj(h1_2),    adj(h1_1_rv)),
  `(3)` = c(adj(h1_3),    adj(h1_1_rv)),
  `(4)` = c(adj(h1_4),    adj(h1_4_rv)),
  `(5)` = c(adj(h1_5),    adj(h1_5_rv))
)
names(test_h1)[2:6] <- c("(1)", "(2)", "(3)", "(4)", "(5)")

stargazer(test_h1, summary = FALSE, type = "html", rownames = FALSE,
          out = 'Table A.1.html')

# --- Table 1 ---
models <- list("(1)" = h1_1, "(2)" = h1_2, "(3)" = h1_3, "(4)" = h1_4,
               "(5)" = h1_5, "(6)" = h1_6, "(7)" = h1_7, "(8)" = h1_16,
               "(9)" = h1_18, "(10)" = h1_19)
modelsummary(models, stars = TRUE, escape = FALSE,
             output = 'Table 1.html',
             coef_map = cm_h1, fmt = 2, vcov = "robust",
             gof_omit = "Std.Errors|AIC|BIC|RMSE|R2 Within|R2 Within Adjusted",
             notes = list("Robust standard errors are clustered by country."))

# --- Table A.2 ---
models <- list("(1)" = h1_8, "(2)" = h1_9, "(3)" = h1_10, "(4)" = h1_11, "(5)" = h1_12)
modelsummary(models, stars = TRUE, escape = FALSE,
             output = 'Table A.2.html',
             coef_map = cm_h1, fmt = 2, vcov = "robust",
             gof_omit = "Std.Errors|AIC|BIC|RMSE|R2 Within|R2 Within Adjusted",
             notes = list("Robust standard errors are clustered by country."))

# --- Table A.3 ---
models <- list("(1)" = h1_n1, "(2)" = h1_n2, "(3)" = h1_n3,
               "(4)" = h1_n4, "(5)" = h1_n5, "(6)" = h1_n6)
modelsummary(models, stars = TRUE, escape = FALSE,
             output = 'Table A.3.html',
             coef_map = cm_h1, fmt = 2, vcov = "robust",
             gof_omit = "Std.Errors|AIC|BIC|RMSE|R2 Within|R2 Within Adjusted",
             notes = list("Robust standard errors are clustered by country."))

################################################################################
# FIGURES
################################################################################

# --- Figure A.1: Correlation matrix ---
corr_vars <- df[, c("BIOs", "Prison", "Police", "Intelligence", "Military",
                    "Paramilitary", "Torture (CIRI)", "Torture (V-Dem)")]

png('Figure A.1.png',
    3000, 3000, res = 475)
ggcorrplot(round(cor(na.omit(corr_vars)), 2),
           method = "square", type = "lower", outline.col = "gray",
           ggtheme = ggplot2::theme_gray,
           colors = c("blue", "white", "red"),
           lab = TRUE, lab_size = 4) +
  ggtitle("") + theme_bw() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, size = 10, hjust = 1),
        axis.text.y = element_text(angle = 0, vjust = 1, size = 10, hjust = 1)) +
  xlab("") + ylab("")
dev.off()

# --- Figure A.4: Distribution of BIO torture observations ---
png('Figure A.4.png',
    3000, 2000, res = 500)
ggplot(df[df$amelia_fit > 0, ], aes(x = amelia_fit)) +
  geom_bar(color = "steelblue", fill = "white") +
  labs(x = "Level of Torture", y = "Frequency") +
  theme_minimal() +
  theme(axis.text.x = element_text(hjust = 1),
        panel.background = element_rect(fill = "white"))
dev.off()

################################################################################
# MAIN OLS MODELS (Table 2 + robustness variants)
################################################################################

# --- Wall construction (Table 2) ---
m1  <- feols(log(amelia_fit+1) ~ wall_const | ccode + year, df, cluster = "ccode")
m2  <- feols(log(amelia_fit+1) ~ wall_const + mean_anxiety + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m3  <- feols(log(amelia_fit+1) ~ wall_const + mean_anger + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m4  <- feols(log(amelia_fit+1) ~ wall_const*mean_anxiety + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m5  <- feols(log(amelia_fit+1) ~ wall_const*mean_anger + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m6  <- feols(log(amelia_fit+1) ~ wall_const + log(migration_intpl+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m7  <- feols(log(amelia_fit+1) ~ wall_const + log(refugees+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m8  <- feols(log(amelia_fit+1) ~ wall_const + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m9  <- feols(log(amelia_fit+1) ~ wall_const*mean_anxiety + log(migration_intpl+1) + log(refugees+1) + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m10 <- feols(log(amelia_fit+1) ~ wall_const*mean_anger + log(migration_intpl+1) + log(refugees+1) + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")

# --- Table A.5: Asylum seeker robustness ---
m7.a  <- feols(log(amelia_fit+1) ~ wall_const + log(asylum_seekers+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m9.a  <- feols(log(amelia_fit+1) ~ wall_const*mean_anxiety + log(asylum_seekers+1) + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m10.a <- feols(log(amelia_fit+1) ~ wall_const*mean_anger + log(asylum_seekers+1) + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")

# --- Table A.6: Border infrastructure ---
m21 <- feols(log(amelia_fit+1) ~ structure_tot | ccode + year, df, cluster = "ccode")
m22 <- feols(log(amelia_fit+1) ~ structure_tot + mean_anxiety + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m23 <- feols(log(amelia_fit+1) ~ structure_tot + mean_anger + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m24 <- feols(log(amelia_fit+1) ~ structure_tot*mean_anxiety + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m25 <- feols(log(amelia_fit+1) ~ structure_tot*mean_anger + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m26 <- feols(log(amelia_fit+1) ~ structure_tot + log(migration_intpl+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m27 <- feols(log(amelia_fit+1) ~ structure_tot + log(refugees+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m28 <- feols(log(amelia_fit+1) ~ structure_tot + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m29 <- feols(log(amelia_fit+1) ~ structure_tot*mean_anxiety + log(migration_intpl+1) + log(refugees+1) + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m30 <- feols(log(amelia_fit+1) ~ structure_tot*mean_anger + log(migration_intpl+1) + log(refugees+1) + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")

# --- Table A.7: Visa extended ---
m41 <- feols(log(amelia_fit+1) ~ visa_extended | ccode + year, df, cluster = "ccode")
m42 <- feols(log(amelia_fit+1) ~ visa_extended + mean_anxiety + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m43 <- feols(log(amelia_fit+1) ~ visa_extended + mean_anger + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m44 <- feols(log(amelia_fit+1) ~ visa_extended*mean_anxiety + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m45 <- feols(log(amelia_fit+1) ~ visa_extended*mean_anger + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m46 <- feols(log(amelia_fit+1) ~ visa_extended + log(migration_intpl+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m47 <- feols(log(amelia_fit+1) ~ visa_extended + log(refugees+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m48 <- feols(log(amelia_fit+1) ~ visa_extended + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m49 <- feols(log(amelia_fit+1) ~ visa_extended*mean_anxiety + log(migration_intpl+1) + log(refugees+1) + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m50 <- feols(log(amelia_fit+1) ~ visa_extended*mean_anger + log(migration_intpl+1) + log(refugees+1) + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")

# --- Table A.8: Wall + border infrastructure ---
bo1 <- feols(log(amelia_fit+1) ~ wall_const + structure_tot | ccode + year, df, cluster = "ccode")
bo2 <- feols(log(amelia_fit+1) ~ wall_const + structure_tot + torture_vdem + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
bo3 <- feols(log(amelia_fit+1) ~ wall_const + structure_tot + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
bo4 <- feols(log(amelia_fit+1) ~ wall_const + structure_tot + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")

# --- Table A.9: Wall x CSO repression ---
visa5 <- feols(log(amelia_fit+1) ~ wall_const*v2csreprss | ccode + year, df, cluster = "ccode")
visa6 <- feols(log(amelia_fit+1) ~ wall_const*v2csreprss + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
visa7 <- feols(log(amelia_fit+1) ~ wall_const*v2csreprss + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
visa8 <- feols(log(amelia_fit+1) ~ wall_const*v2csreprss + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")

# --- Table 2 output ---
models <- list("(1)" = m1, "(2)" = m2, "(3)" = m3, "(4)" = m4, "(5)" = m5,
               "(6)" = m6, "(7)" = m7, "(8)" = m8, "(9)" = m9, "(10)" = m10)
modelsummary(models, stars = TRUE, escape = FALSE,
             output = 'Table 2.html',
             coef_map = cm_main, fmt = 2, vcov = "robust",
             gof_omit = "Std.Errors|FE|AIC|BIC|RMSE|R2 Within|R2 Within Adjusted",
             notes = list("Robust standard errors are clustered by country. Country and Year FE included."))

# --- Table A.5 ---
models <- list("(1)" = m7.a, "(2)" = m9.a, "(3)" = m10.a)
modelsummary(models, stars = TRUE, escape = FALSE,
             output = "Table A.5.html",
             coef_map = cm_main, fmt = 2, vcov = "robust",
             gof_omit = "Std.Errors|FE|AIC|BIC|RMSE|R2 Within|R2 Within Adjusted",
             notes = list("Robust standard errors are clustered by country. Country and Year FE included."))

# --- Table A.6 ---
models <- list("(1)" = m21, "(2)" = m22, "(3)" = m23, "(4)" = m24, "(5)" = m25,
               "(6)" = m26, "(7)" = m27, "(8)" = m28, "(9)" = m29, "(10)" = m30)
modelsummary(models, stars = TRUE, escape = FALSE,
             output = "Table A.6.html",
             coef_map = cm_main, fmt = 2, vcov = "robust",
             gof_omit = "Std.Errors|FE|AIC|BIC|RMSE|R2 Within|R2 Within Adjusted",
             notes = list("Robust standard errors are clustered by country. Country and Year FE included."))

# --- Table A.7 ---
models <- list("(1)" = m41, "(2)" = m42, "(3)" = m43, "(4)" = m44, "(5)" = m45,
               "(6)" = m46, "(7)" = m47, "(8)" = m48, "(9)" = m49, "(10)" = m50)
modelsummary(models, stars = TRUE, escape = FALSE,
             output = "Table A.7.html",
             coef_map = cm_main, fmt = 2, vcov = "robust",
             gof_omit = "Std.Errors|FE|AIC|BIC|RMSE|R2 Within|R2 Within Adjusted",
             notes = list("Robust standard errors are clustered by country. Country and Year FE included."))

# --- Table A.8 ---
models <- list("(1)" = bo1, "(2)" = bo2, "(3)" = bo3, "(4)" = bo4)
modelsummary(models, stars = TRUE, escape = FALSE,
             output = 'Table A.8.html',
             coef_map = cm_main, fmt = 3, vcov = "robust",
             gof_omit = "Std.Errors|FE|AIC|BIC|RMSE|R2 Within|R2 Within Adjusted",
             notes = list("Robust standard errors are clustered by country. Country and Year FE included."))

# --- Table A.9 ---
models <- list("(1)" = visa5, "(2)" = visa6, "(3)" = visa7, "(4)" = visa8)
modelsummary(models, stars = TRUE, escape = FALSE,
             output = "Table A.9.html",
             coef_map = cm_main, fmt = 3, vcov = "robust",
             gof_omit = "Std.Errors|FE|AIC|BIC|RMSE|R2 Within|R2 Within Adjusted",
             notes = list("Robust standard errors are clustered by country. Country and Year FE included."))

################################################################################
# PLACEBO TESTS (Table A.14)
################################################################################

pl1 <- feols(Police ~ wall_const + torture_vdem + log(wall_reporting+1) + ideology + civil_war + log(fatalities+1) + liberal_democracy + log(population) + log(gdp_pc) + log(migration_intpl+1) | ccode + year, data = df, cluster = "ccode")
pl2 <- feols(Prison ~ wall_const + torture_vdem + log(wall_reporting+1) + ideology + civil_war + log(fatalities+1) + liberal_democracy + log(population) + log(gdp_pc) + log(migration_intpl+1) | ccode + year, data = df, cluster = "ccode")
pl3 <- feols(Military ~ wall_const + torture_vdem + log(wall_reporting+1) + ideology + civil_war + log(fatalities+1) + liberal_democracy + log(population) + log(gdp_pc) + log(migration_intpl+1) | ccode + year, data = df, cluster = "ccode")
pl4 <- feols(Intelligence ~ wall_const + torture_vdem + log(wall_reporting+1) + ideology + civil_war + log(fatalities+1) + liberal_democracy + log(population) + log(gdp_pc) + log(migration_intpl+1) | ccode + year, data = df, cluster = "ccode")
pl5 <- feols(Paramilitary ~ wall_const + torture_vdem + log(wall_reporting+1) + ideology + civil_war + log(fatalities+1) + liberal_democracy + log(population) + log(gdp_pc) + log(migration_intpl+1) | ccode + year, data = df, cluster = "ccode")

models <- list("Police" = pl1, "Prison" = pl2, "Military" = pl3, "Intelligence" = pl4, "Paramilitary" = pl5)
modelsummary(models, stars = TRUE, escape = FALSE,
             output = 'Table A.14.html',
             coef_map = cm_main, fmt = 2, vcov = "robust",
             gof_omit = "Std.Errors|AIC|BIC|RMSE|R2 Within|R2 Within Adjusted",
             notes = list("Robust standard errors are clustered by country."))

################################################################################
# ROBUSTNESS — LOGIT (Table A.10)
################################################################################

l1  <- feglm(torture_dum ~ wall_const | ccode + year, df, cluster = "ccode", family = binomial(link = "logit"))
l2  <- feglm(torture_dum ~ wall_const + mean_anxiety + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = binomial(link = "logit"))
l3  <- feglm(torture_dum ~ wall_const + mean_anger + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = binomial(link = "logit"))
l4  <- feglm(torture_dum ~ wall_const*mean_anxiety + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = binomial(link = "logit"))
l5  <- feglm(torture_dum ~ wall_const*mean_anger + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = binomial(link = "logit"))
l6  <- feglm(torture_dum ~ wall_const + log(migration_intpl+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = binomial(link = "logit"))
l7  <- feglm(torture_dum ~ wall_const + log(refugees+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = binomial(link = "logit"))
l8  <- feglm(torture_dum ~ wall_const + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = binomial(link = "logit"))
l9  <- feglm(torture_dum ~ wall_const*mean_anxiety + log(migration_intpl+1) + log(refugees+1) + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = binomial(link = "logit"))
l10 <- feglm(torture_dum ~ wall_const*mean_anger + log(migration_intpl+1) + log(refugees+1) + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = binomial(link = "logit"))

models <- list("(1)" = l1, "(2)" = l2, "(3)" = l3, "(4)" = l4, "(5)" = l5,
               "(6)" = l6, "(7)" = l7, "(8)" = l8, "(9)" = l9, "(10)" = l10)
modelsummary(models, stars = TRUE, escape = FALSE,
             output = 'Table A.10.html',
             coef_map = cm_main, fmt = 2, vcov = "robust",
             gof_omit = "Std.Errors|AIC|BIC|RMSE|R2 Within|R2 Within Adjusted",
             notes = list("Robust standard errors are clustered by country."))

################################################################################
# ROBUSTNESS — POISSON (Table A.11) & NEGATIVE BINOMIAL (Table A.12)
################################################################################

pois1  <- feglm(amelia_fit ~ wall_const | ccode + year, df, cluster = "ccode", family = "poisson")
pois2  <- feglm(amelia_fit ~ wall_const + mean_anxiety + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = "poisson")
pois3  <- feglm(amelia_fit ~ wall_const + mean_anger + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = "poisson")
pois4  <- feglm(amelia_fit ~ wall_const*mean_anxiety + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = "poisson")
pois5  <- feglm(amelia_fit ~ wall_const*mean_anger + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = "poisson")
pois6  <- feglm(amelia_fit ~ wall_const + log(migration_intpl+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = "poisson")
pois7  <- feglm(amelia_fit ~ wall_const + log(refugees+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = "poisson")
pois8  <- feglm(amelia_fit ~ wall_const + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = "poisson")
pois9  <- feglm(amelia_fit ~ wall_const*mean_anxiety + log(migration_intpl+1) + log(refugees+1) + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = "poisson")
pois10 <- feglm(amelia_fit ~ wall_const*mean_anger + log(migration_intpl+1) + log(refugees+1) + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode", family = "poisson")

models <- list("(1)" = pois1, "(2)" = pois2, "(3)" = pois3, "(4)" = pois4, "(5)" = pois5,
               "(6)" = pois6, "(7)" = pois7, "(8)" = pois8, "(9)" = pois9, "(10)" = pois10)
modelsummary(models, stars = TRUE, escape = FALSE,
             output = 'Table A.11.html',
             coef_map = cm_main, fmt = 2, vcov = "robust",
             gof_omit = "Std.Errors|AIC|BIC|RMSE|R2.Adj|R2 Within|R2 Within Adjusted",
             notes = list("Robust standard errors are clustered by country."))

negbi1  <- fixest::femlm(amelia_fit ~ wall_const | ccode + year, df, family = "negbin", panel.id = c("ccode", "year"))
negbi2  <- fixest::femlm(amelia_fit ~ wall_const + mean_anxiety + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, family = "negbin", panel.id = c("ccode", "year"))
negbi3  <- fixest::femlm(amelia_fit ~ wall_const + mean_anger + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, family = "negbin", panel.id = c("ccode", "year"))
negbi4  <- fixest::femlm(amelia_fit ~ wall_const*mean_anxiety + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, family = "negbin", panel.id = c("ccode", "year"))
negbi5  <- fixest::femlm(amelia_fit ~ wall_const*mean_anger + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, family = "negbin", panel.id = c("ccode", "year"))
negbi6  <- fixest::femlm(amelia_fit ~ wall_const + log(migration_intpl+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, family = "negbin", panel.id = c("ccode", "year"))
negbi7  <- fixest::femlm(amelia_fit ~ wall_const + log(refugees+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, family = "negbin", panel.id = c("ccode", "year"))
negbi8  <- fixest::femlm(amelia_fit ~ wall_const + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, family = "negbin", panel.id = c("ccode", "year"))
negbi9  <- fixest::femlm(amelia_fit ~ wall_const*mean_anxiety + log(migration_intpl+1) + log(refugees+1) + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, family = "negbin", panel.id = c("ccode", "year"))
negbi10 <- fixest::femlm(amelia_fit ~ wall_const*mean_anger + log(migration_intpl+1) + log(refugees+1) + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, family = "negbin", panel.id = c("ccode", "year"))

models <- list("(1)" = negbi1, "(2)" = negbi2, "(3)" = negbi3, "(4)" = negbi4, "(5)" = negbi5,
               "(6)" = negbi6, "(7)" = negbi7, "(8)" = negbi8, "(9)" = negbi9, "(10)" = negbi10)
modelsummary(models, stars = TRUE, escape = FALSE,
             output = "Table A.12.html",
             coef_map = cm_main, fmt = 2, vcov = "robust",
             gof_omit = "Std.Errors|AIC|BIC|RMSE|R2.Adj.|R2 Within|R2 Within Adjusted",
             notes = list("Robust standard errors are clustered by country."))

################################################################################
# ROBUSTNESS — 3-Year EFFECT (Table A.13)
################################################################################

m1_3  <- feols(log(amelia_fit+1) ~ wall_const3 | ccode + year, df, cluster = "ccode")
m2_3  <- feols(log(amelia_fit+1) ~ wall_const3 + mean_anxiety + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m3_3  <- feols(log(amelia_fit+1) ~ wall_const3 + mean_anger + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m4_3  <- feols(log(amelia_fit+1) ~ wall_const3*mean_anxiety + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m5_3  <- feols(log(amelia_fit+1) ~ wall_const3*mean_anger + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m6_3  <- feols(log(amelia_fit+1) ~ wall_const3 + log(migration_intpl+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m7_3  <- feols(log(amelia_fit+1) ~ wall_const3 + log(refugees+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m8_3  <- feols(log(amelia_fit+1) ~ wall_const3 + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m9_3  <- feols(log(amelia_fit+1) ~ wall_const3*mean_anxiety + log(migration_intpl+1) + log(refugees+1) + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")
m10_3 <- feols(log(amelia_fit+1) ~ wall_const3*mean_anger + log(migration_intpl+1) + log(refugees+1) + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year, df, cluster = "ccode")

models <- list("(1)" = m1_3, "(2)" = m2_3, "(3)" = m3_3, "(4)" = m4_3, "(5)" = m5_3,
               "(6)" = m6_3, "(7)" = m7_3, "(8)" = m8_3, "(9)" = m9_3, "(10)" = m10_3)
modelsummary(models, stars = TRUE, escape = FALSE,
             output = "Table A.13.html",
             coef_map = cm_main, fmt = 2, vcov = "robust",
             gof_omit = "Std.Errors|AIC|BIC|RMSE|R2 Within|R2 Within Adjusted",
             notes = list("Robust standard errors are clustered by country."))

################################################################################
# JACKNIFE (Figure A.2)
################################################################################

ccode_s      <- unique(df$ccode)
observations2 <- length(ccode_s)
jack.reg2    <- numeric(observations2)
se2          <- numeric(observations2)

for (i in 1:observations2) {
  model <- feols(log(amelia_fit+1) ~ wall_const + torture_vdem + ideology + log(wall_reporting+1) + log(migration_intpl+1) + log(fatalities+1) + liberal_democracy + log(population) | ccode + year,
                 df[df$ccode != ccode_s[i], ], cluster = "ccode")
  jack.reg2[i] <- coef(model)[1]
  se2[i]       <- summary(model)$se[1]
  print(paste0(round((i / observations2) * 100, 2), "% completed."))
}

jacknife      <- data.frame(country_out = ccode_s, est = jack.reg2, se = se2)
jacknife$lb   <- jacknife$est - (1.96  * jacknife$se)
jacknife$ub   <- jacknife$est + (1.96  * jacknife$se)
jacknife$lb2  <- jacknife$est - (1.645 * jacknife$se)
jacknife$ub2  <- jacknife$est + (1.645 * jacknife$se)
jacknife$group <- ifelse(jacknife$country_out < 451, "", " ")

png('Figure A.2.png',
    3000, 2500, res = 350)
ggplot(jacknife, aes(x = est, y = as.factor(country_out))) +
  geom_point(position = position_dodge(width = 0.5)) +
  geom_errorbarh(aes(xmin = lb, xmax = ub),   height = 0, size = 0.35, position = position_dodge(width = 0.5)) +
  geom_errorbarh(aes(xmin = lb2, xmax = ub2), height = 0, size = 0.65, position = position_dodge(width = 0.5)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
  labs(x = "Wall Coefficient Estimate", y = "Country Out") +
  theme_minimal() +
  facet_wrap(~ group, ncol = 2, scales = "free_y")
dev.off()

################################################################################
# SENSITIVITY ANALYSIS
################################################################################

band <- c(1, 20, 50)

sens2a <- sensemakr::sensemakr(model = m2, treatment = "wall_const", benchmark_covariates = "torture_vdem",    kd = band, q = 5, alpha = 0.95, reduce = TRUE)
sens2b <- sensemakr::sensemakr(model = m2, treatment = "wall_const", benchmark_covariates = "liberal_democracy", kd = band, q = 5, alpha = 0.95, reduce = TRUE)
sens2c <- sensemakr::sensemakr(model = m2, treatment = "wall_const", benchmark_covariates = "log(population)", kd = band, q = 5, alpha = 0.95, reduce = TRUE)
sens2d <- sensemakr::sensemakr(model = m2, treatment = "wall_const", benchmark_covariates = "ideology",        kd = band, q = 5, alpha = 0.95, reduce = TRUE)
sens2e <- sensemakr::sensemakr(model = m2, treatment = "wall_const", benchmark_covariates = "wall_reporting",  kd = band, q = 5, alpha = 0.95, reduce = TRUE)
sens6  <- sensemakr::sensemakr(model = m6, treatment = "wall_const", benchmark_covariates = "log(migration_intpl + 1)", kd = band, q = 5, alpha = 0.95, reduce = TRUE)
sens7  <- sensemakr::sensemakr(model = m7, treatment = "wall_const", benchmark_covariates = "log(refugees + 1)",        kd = band, q = 5, alpha = 0.95, reduce = TRUE)
sens8a <- sensemakr::sensemakr(model = m8, treatment = "wall_const", benchmark_covariates = "civil_war",       kd = band, q = 5, alpha = 0.95, reduce = TRUE)
sens8b <- sensemakr::sensemakr(model = m8, treatment = "wall_const", benchmark_covariates = "log(fatalities + 1)", kd = band, q = 5, alpha = 0.95, reduce = TRUE)

sens2a$bounds[, 1] <- c("1x Torture",       "20x Torture",       "50x Torture")
sens2b$bounds[, 1] <- c("1x Lib. Dem",      "20x Lib. Dem",      "50x Lib. Dem")
sens2c$bounds[, 1] <- c("1x Population",    "20x Population",    "50x Population")
sens2d$bounds[, 1] <- c("1x Ideology",      "20x Ideology",      "50x Ideology")
sens2e$bounds[, 1] <- c("1x Wall Reporting","20x Wall Reporting","50x Wall Reporting")
sens6$bounds[, 1]  <- c("1x Migrants",      "20x Migrants",      "50x Migrants")
sens7$bounds[, 1]  <- c("1x Refugees",      "20x Refugees",      "50x Refugees")
sens8a$bounds[, 1] <- c("1x Civil War",     "20x Civil War",     "50x Civil War")
sens8b$bounds[, 1] <- c("1x Violence",      "20x Violence",      "50x Violence")

sens_list <- list(
  list(obj = sens2a, file = "1"), list(obj = sens2b, file = "2"),
  list(obj = sens2c, file = "3"), list(obj = sens2d, file = "4"),
  list(obj = sens2e, file = "5"), list(obj = sens6,  file = "6"),
  list(obj = sens7,  file = "7"),  list(obj = sens8a, file = "8"),
  list(obj = sens8b, file = "9")
)

################################################################################
# FIGURE A.3
################################################################################

for (s in sens_list) {
  png(paste0('Figure A.3_', s$file, '.png'),
      2400, 2400, res = 500)
  plot(s$obj, sensitivity.of = "t-value", cex.label.text = 0.5)
  dev.off()
}

################################################################################
# INSTRUMENTAL VARIABLE (Tables 3 & A.16)
################################################################################

invt_min_1 <- abs(min(na.omit(df$nonfinancial))) + 1

iv1  <- feols(log(amelia_fit+1) ~ 1 | ccode + year | wall_const ~ log(nonfinancial+invt_min_1), df, cluster = "ccode")
iv2  <- feols(log(amelia_fit+1) ~ mean_anxiety + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year | wall_const ~ log(nonfinancial+invt_min_1), df, cluster = "ccode")
iv3  <- feols(log(amelia_fit+1) ~ mean_anger + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year | wall_const ~ log(nonfinancial+invt_min_1), df, cluster = "ccode")
iv4  <- feols(log(amelia_fit+1) ~ log(migration_intpl+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year | wall_const ~ log(nonfinancial+invt_min_1), df, cluster = "ccode")
iv5  <- feols(log(amelia_fit+1) ~ log(refugees+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year | wall_const ~ log(nonfinancial+invt_min_1), df, cluster = "ccode")
iv6  <- feols(log(amelia_fit+1) ~ civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year | wall_const ~ log(nonfinancial+invt_min_1), df, cluster = "ccode")
iv7 <- feols(log(amelia_fit+1) ~ mean_anxiety + log(refugees+1) + log(migration_intpl+1) + civil_war + log(fatalities+1) + ideology + torture_vdem + liberal_democracy + wall_reporting + log(population) | ccode + year | wall_const ~ log(nonfinancial+invt_min_1), df, cluster = "ccode")

# Table 3
models <- list("(1)" = iv1, "(2)" = iv2, "(3)" = iv3, "(4)" = iv4,
               "(5)" = iv5, "(6)" = iv6, "(7)" = iv7)
modelsummary(models, stars = TRUE, escape = FALSE,
             output = "Table 3.html",
             coef_map = cm_main, fmt = 2, vcov = "robust",
             gof_omit = "Std.Errors|AIC|BIC|RMSE|R2 Within|R2 Within Adjusted",
             notes = list("Robust standard errors included. Country and Year FE included."))

# Anderson-Rubin test p-values via ivDiag with no bootstrapping (near-instant)
ivdiag1 <- ivDiag(df, Y="amelia_log", D="wall_const", Z="nonfinancial_log", controls = NULL,c("ccode","year"), nboots = 1000, seed = 20240512)
closeAllConnections()
ivdiag2 <- ivDiag(df, Y="amelia_log", D="wall_const", Z="nonfinancial_log", controls = c("mean_anxiety","ideology","torture_vdem","liberal_democracy","wall_log","population_log"), FE = c("ccode","year"), cl = "ccode", nboots = 1000, seed = 20240512)
closeAllConnections()
ivdiag3 <- ivDiag(df, Y="amelia_log", D="wall_const", Z="nonfinancial_log", controls = c("mean_anger","ideology","torture_vdem","liberal_democracy","wall_log","population_log"), FE = c("ccode","year"), cl = "ccode", nboots = 1000, seed = 20240512)
closeAllConnections()
ivdiag4 <- ivDiag(df, Y="amelia_log", D="wall_const", Z="nonfinancial_log", controls = c("migrants_log","ideology","torture_vdem","liberal_democracy","wall_log","population_log"), FE = c("ccode","year"), cl = "ccode", nboots = 1000, seed = 20240512)
closeAllConnections()
ivdiag5 <- ivDiag(df, Y="amelia_log", D="wall_const", Z="nonfinancial_log", controls = c("refugees_log","ideology","torture_vdem","liberal_democracy","wall_log","population_log"), FE = c("ccode","year"), cl = "ccode", nboots = 1000, seed = 20240512)
closeAllConnections()
ivdiag6 <- ivDiag(df, Y="amelia_log", D="wall_const", Z="nonfinancial_log", controls = c("civil_war","fatalities_log","ideology","torture_vdem","liberal_democracy","wall_log","population_log"), FE = c("ccode","year"), cl = "ccode", nboots = 1000, seed = 20240512)
closeAllConnections()
ivdiag7 <- ivDiag(df, Y="amelia_log", D="wall_const", Z="nonfinancial_log", controls = c("mean_anxiety","migrants_log","refugees_log","civil_war","fatalities_log","ideology","torture_vdem","liberal_democracy","wall_log","population_log"), FE = c("ccode","year"), cl = "ccode", nboots = 1000, seed = 20240512)

# Table A.16 (first stage + diagnostics)
row <- data.frame("Coefficients" = c("F-statistic","IV Wald test","Anderson-Rubin test"),
                  "(1)"=c(fitstat(iv1,"ivf")$`ivf1::wall_const`$stat,fitstat(iv1,"ivwald")$`ivwald1::wall_const`$p,ivdiag1$AR$Fstat[4]),
                  "(2)"=c(fitstat(iv2,"ivf")$`ivf1::wall_const`$stat,fitstat(iv2,"ivwald")$`ivwald1::wall_const`$p,ivdiag2$AR$Fstat[4]),
                  "(3)"=c(fitstat(iv3,"ivf")$`ivf1::wall_const`$stat,fitstat(iv3,"ivwald")$`ivwald1::wall_const`$p,ivdiag3$AR$Fstat[4]),
                  "(4)"=c(fitstat(iv4,"ivf")$`ivf1::wall_const`$stat,fitstat(iv4,"ivwald")$`ivwald1::wall_const`$p,ivdiag4$AR$Fstat[4]),
                  "(5)"=c(fitstat(iv5,"ivf")$`ivf1::wall_const`$stat,fitstat(iv5,"ivwald")$`ivwald1::wall_const`$p,ivdiag5$AR$Fstat[4]),
                  "(6)"=c(fitstat(iv6,"ivf")$`ivf1::wall_const`$stat,fitstat(iv6,"ivwald")$`ivwald1::wall_const`$p,ivdiag6$AR$Fstat[4]),
                  "(7)"=c(fitstat(iv7,"ivf")$`ivf1::wall_const`$stat,fitstat(iv7,"ivwald")$`ivwald1::wall_const`$p,ivdiag7$AR$Fstat[4]))

attr(row, "position") <- c(27,28,29)

models <- list("(1)" = summary(iv1,  stage = 1), "(2)" = summary(iv2,  stage = 1),
               "(3)" = summary(iv3,  stage = 1), "(4)" = summary(iv4,  stage = 1),
               "(5)" = summary(iv5,  stage = 1), "(6)" = summary(iv6,  stage = 1),
               "(7)" = summary(iv7, stage = 1))
modelsummary(models, add_rows = row, stars = TRUE, escape = FALSE,
             output = "Table A.16.html",
             coef_map = cm_main, fmt = 2, vcov = "robust",
             gof_omit = "Std.Errors|AIC|BIC|RMSE|R2 Within|R2 Within Adjusted|R2.Adj",
             notes = list("Robust standard errors included. Country and Year FE included."))

################################################################################
# DIFFERENCE-IN-DIFFERENCES (Table A.15 & Figure 3)
################################################################################

mw.attgt <- att_gt(
  yname      = "amelia_log",
  gname      = "wall_construction2",
  idname     = "ccode",
  tname      = "year",
  data       = df,
  allow_unbalanced_panel = TRUE,
  control_group = "notyettreated",
  biters     = 10000,
  base_period = "universal",
  pl         = TRUE,
  clustervars = "ccode"
)

agg.es <- aggte(mw.attgt, type = "dynamic", na.rm = TRUE, alp = 0.95,
                min_e = -7, max_e = 7, clustervars = "ccode")

# Table A.15
# Stars need to be added manually.
modelsummary(list("(1)" = agg.es), stars = TRUE, escape = FALSE,
             output = 'Table A.15.html',
             vcov = "classical", fmt = 2,
             gof_omit = "Std.Errors|AIC|BIC|RMSE")

# Figure 3
png('Figure 3.png',
    2500, 1700, res = 400)
ggdid(agg.es) +
  xlab("Time since Wall Construction") +
  ylab("Estimated Level of Torture change") +
  ggtitle("") +
  theme(axis.text = element_text(size = 5))
dev.off()